/*
 *  lattice.cpp
 *  EPPI-0.0
 *
 *  Created by Juan Carlos Garcia Orden on 02/02/09.
 *
 *  modified by Markus Lutz on 26/10/09
 */

#include "lattice.h"

Lattice::Lattice(ifstream &in, string cname, MecSystem *pSys)
:	Body(in,cname,pSys),
	meshFileName(),
	nNodes(),
	nElements(),
//	vPointNames(),
//	connect(),
	connectSurface(), // edited
	vCoord(),
	vertexTag(),
	nVertex(),
	matName(""),
	area(1.0)
{	
	dim = parentSystem->dim;	/** Must be set always.  dim belongs to the mother class Body */
	type = "lattice";
	
	vCoord.resize(dim);
	vector< vector<string> > connect;
	
 	ifstream ini;
	char keyword[20];
	string auxs;				/** Aux string */
	vector<string> auxsv;		/** Aux string vector */
	vector<string> vertexNames;
	char contactFileName[30];
	int numContact;
	string line;
	vector<string> vPNames;
	int i;
	ifstream inContact;	

	while (in >> keyword && strcmp(keyword,"ENDBODY") )
	{
		if (!strcmp(keyword,"MESHFILE"))
		{
			/** keyword has MESHFILE */
			in >> meshFileName;	/** Reads the name of the file */
		
			in >> auxs;				/** Reads the AREA label */
			in >> area;				/** Reads the area value */
			in >> auxs;				/** Reads the MAT label */
			in >> matName;			/** Reads the material name */

			ini.open(meshFileName);
			// gmsh format
			readGmsh(ini,connect);
			getBoundary(connect); // edit
			ini.close();
		}
		if (!strcmp(keyword,"CONTACTPOINTS"))
		{
			in >> auxs;
			if (!auxs.compare("all"))	contactAllOption = true;
			if (!auxs.compare("file"))
			{
				in >> contactFileName;
				contactFileOption = true;
			}
			if (!auxs.compare("names"))
			{
				contactNamesOption = true;
				getline(in,line);
				istringstream iss(line);
				while (iss >> keyword )
				{
					auxs = keyword;
					vPNames.push_back(auxs);
					//cout << "keyword is " << keyword << endl;
				}				
			}
		}
	}
	if(parentSystem->getProcNum()==parentSystem->getRankNum() ) cout << "	Body " << name << " succesfully read \n";
		
	// If nothing is said about contact, take by default the "all" option
	if (contactAllOption==false && contactFileOption==false && contactNamesOption==false)	contactAllOption = true;
	
	if (contactAllOption==true)
	{
		vPointNamesContact = vPointNamesSurface;
	}
	if (contactFileOption==true)
	{
		inContact.open(contactFileName);
		inContact >> numContact;
		for (i=0; i<numContact; ++i)
		{
			inContact >> auxs;
			vPointNamesContact.push_back(auxs);
		}
		inContact.close();
	}
	if (contactNamesOption==true)
	{
		for (i=0; i<vPNames.size(); ++i)	vPointNamesContact.push_back(vPNames[i]);
	}
}

// ***

void Lattice::setPTPointsComp()
{
	int j, k;
	vector<string> auxs;
	map<string, Point*, less<string> > maux;
	pair<string,Point*> mauxEntry;
	pair<string,int> sharedEntry;
	
	for(j=0; j<vComponents.size(); j++)		
	//Loop on the components
	{
		auxs = (vComponents[j]->getPointNames());
		for(k=0; k<auxs.size(); k++)		
		// Loop on the names of the points within each component
		{
			/** If the point it is not at the local process ... */
			if( (parentSystem->checkIfPointProcess(auxs[k]))==false )
			{
				/** ...  create it at root system, copying the original point from the remote process */
				parentSystem->createPointAtRoot( parentSystem->getPTPoint(auxs[k]) );
				/** ...  and add the name of the shared point at the remote process */
				parentSystem->addRemoteSharedPoint(auxs[k],parentSystem->getProcNum());
			}
			mauxEntry.first = auxs[k];
			mauxEntry.second = (parentSystem->getPTPointLocal(auxs[k]));
			maux.insert(mauxEntry);
		}
		vComponents[j]->setMPoints(maux);
		maux.clear();
	}
}

// ***

void Lattice::setPTPointsLoad()
{
	int j, k;
	vector<string> auxs;
	map<string, Point*, less<string> > maux;
	pair<string,Point*> mauxEntry;
	pair<string,int> sharedEntry;
	
	for(j=0; j<vLoads.size(); j++)		
	//Loop on the Loads
	{
		auxs = (vLoads[j]->getPointNames());
		for(k=0; k<auxs.size(); k++)		
		// Loop on the names of the points within each component
		{
			/** If the point it is not at the local process ... */
			if( (parentSystem->checkIfPointProcess(auxs[k]))==false )
			{
				/** ...  create it at root system, copying the original point from the remote process */
				parentSystem->createPointAtRoot( parentSystem->getPTPoint(auxs[k]) );
				/** ...  and add the name of the shared point at the remote process */
				parentSystem->addRemoteSharedPoint(auxs[k],parentSystem->getProcNum());
			}
			mauxEntry.first = auxs[k];
			mauxEntry.second = (parentSystem->getPTPointLocal(auxs[k]));
			maux.insert(mauxEntry);
		}
		vLoads[j]->setMPoints(maux);
		maux.clear();
	}
}

// ***

void Lattice::setPPointSurface()
{
	// Build the vector of pointers to the surface points
	vPPointSurface.resize(vPointNamesSurface.size());
	
	int i;
	string auxs;
	//vector<Point*>	vPoints(nNodes);	
	for(i=0;i<vPointNamesSurface.size();i++)
	{
		auxs = vPointNamesSurface[i];
		addSystemName(auxs,parentSystem->name);
		vPPointSurface[i] = parentSystem->getPTPointLocal( auxs ); 
	}			
}

// ***

void Lattice::setMatAtComp()
{
	int i;
	string matName;
	Material* mat;
	for(i=0; i<vComponents.size(); i++)
	{
		matName = vComponents[i]->getMatName();
		mat = parentSystem->getPTMaterial(matName);		
		vComponents[i]->setMatAtComp(mat);
	}
}

// ***

void Lattice::getBoundary(vector<vector<string> >& connect)
{
// edit1: starts here ....



//*********************************************
// converting ...
//*********************************************
	
	switch (vertexTag) // Element tag (2=triangle, 3=quadrangle, 4=tetrahedron, 5=hexahedron)
	{
		case 2: // triangle
		{	
			vector <string> surface(nVertex,"0");	// tempory vector for the triangle + a identifier 

			// now the triangle is reconstructed by 3 lines
			for(int i = 0; i < nElements;  i += 1)	// loop over all Elements
			{
	
		
				// first line
				surface[0] = connect[i][0];
				surface[1] = connect[i][1];
				surface[2] = surface[0] + 'x' + surface[1]; // identifier

				connectSurface.push_back(surface);

				// second line
				surface[0] = connect[i][1];
				surface[1] = connect[i][2];
				surface[2] = surface[0] + 'x' + surface[1]; // identifier

				connectSurface.push_back(surface);

				// third line
				surface[0] = connect[i][2];
				surface[1] = connect[i][0];
				surface[2] = surface[0] + 'x' + surface[1]; // identifier

				connectSurface.push_back(surface);
			}
			break;
		}
		
		case 3: {cout << "cuadro not available in lattice\n"; break;} 

		case 4: // tetrahedron
		{			
			vector <string> surface(nVertex,"0");	// tempory vector for the triangle + a identifier 

			// now the tetrahedron is reconstructed by 3 triangles
			for(int i = 0; i < nElements;  i += 1)	// loop over all Elements
			{	
		
				// first surface
				surface[0] = connect[i][0];
				surface[1] = connect[i][3];
				surface[2] = connect[i][2];
				surface[3] = surface[0] + 'a' + surface[1] + 'b' + surface[2]; // identifier

				connectSurface.push_back(surface);

				// second surface
				surface[0] = connect[i][0];
				surface[1] = connect[i][1];
				surface[2] = connect[i][2];
				surface[3] = surface[0] + 'c' + surface[1] + 'd' + surface[2]; // identifier

				connectSurface.push_back(surface);

				// third surface
				surface[0] = connect[i][1];
				surface[1] = connect[i][2];
				surface[2] = connect[i][3];
				surface[3] = surface[0] + 'e' + surface[1] + 'f' + surface[2]; // identifier

				connectSurface.push_back(surface);

				// fourth surface
				surface[0] = connect[i][0];
				surface[1] = connect[i][3];
				surface[2] = connect[i][1];
				surface[3] = surface[0] + 'g' + surface[1] + 'h' + surface[2]; // identifier

				connectSurface.push_back(surface);
		
			}
			break;
		}

		case 5: { cout << "hexadron not available in lattice\n"; break;	}

		default:{cout << "vertexTag " << vertexTag << " unknown or undefined\n"; break;}
	}	

	


//*******************************************************
// reductant surfaces are deleted
//*******************************************************

	// vectors	
	vector< vector<string> > save_vector = connectSurface; 	// vector which saves the original structure of the surfaces
	vector< vector<string> > temp_vector;			// vector for the elements of the surface
	vector<string> match;					// vector for the element identifiers of the surface


	// the nodes in all elements are sorted, last element of the vector is not considered
	for (int i = 0; i<connectSurface.size();i++){
		sort (connectSurface[i].begin(), connectSurface[i].end()-1);	
	}

	
	// now the vectors containing the nodes are sorted
	sort(connectSurface.begin(),connectSurface.end());


	// as the vectors are sorted now, only the neighbour of each vector
	// has to be compared for equality	
  	for (int i = 0; i<connectSurface.size()-1;i++)	// loop over all vectors
	{
		if (!equal( connectSurface[i].begin(), connectSurface[i].end()-1, connectSurface[i+1].begin())){
			match.push_back(connectSurface[i][nVertex-1]);
		}
		else {i += 1;}
	}

	

	// the last surface is a special case...
	if (!equal( connectSurface[connectSurface.size()-1].begin(), connectSurface[connectSurface.size()-1].end()-1, connectSurface[connectSurface.size()-2].begin()))
	{
		match.push_back(connectSurface[connectSurface.size()-1][nVertex-1]);
	}


	sort(match.begin(),match.end());	// arrange vector

	for (int i = 0; i< connectSurface.size(); i++ ) // loop over all surfaces
	{	
		if (binary_search(match.begin(),match.end(),save_vector[i][nVertex-1])){	// elements of the surface are compared to the vector containing all
			temp_vector.push_back(save_vector[i]);				// "new" elements, if found
		}
	}
	connectSurface.resize(temp_vector.size());
	connectSurface = temp_vector;		// vector with new elements of the surface

//*******************************************************
// reductant nodes are deleted
//*******************************************************
		
	for (int i = 0; i<connectSurface.size();i++) 
	{	// make big vector with all nodes
		for (int j = 0; j < nVertex-1; j++)
		{
			vPointNamesSurface.push_back(connectSurface[i][j]);
		}
	}
		
	sort(vPointNamesSurface.begin(),vPointNamesSurface.end()); 	// arrange nodes 	
		
	vector <string>::iterator it;
	it = unique (vPointNamesSurface.begin(), vPointNamesSurface.end()); 	// delete all multiple nodes in vector
  	vPointNamesSurface.resize( it - vPointNamesSurface.begin() ); 		// resize node vector
	// and edit1 ends here
}

// ***

void Lattice::readGmsh(ifstream &in, vector<vector<string> >& connect)
{
	string auxs;
	vector<string> vVertex;		/** Auxiliary vector with the names of the vertexes in one element */
	int auxi1, auxi2, auxi3;
	int i;
	int j;
	double coo;
	vector<string> auxsv;		/** Aux string vector. Contains two vertexes, for defining the trusses */
	auxsv.resize(2);
	vector<string> vPointNames;
		
	in >> auxs;		/** $MeshFormat */
	in >> auxi1 >> auxi2 >> auxi3;
	in >> auxs;		/** $EndMeshFormat */
	in >> auxs;		/** $Nodes */
	in >> nNodes;
	vPointNames.resize(nNodes);

	for (i=0; i<nNodes; ++i)	// Read coordinates
	{
		in >> auxs;
		vPointNames[i] = auxs;				// Name of the node
		addSystemName(auxs,parentSystem->name);			/** System name is added to the point name */

		in >> coo;			// First coordinate
		vCoord[0] = coo;
		in >> coo;			// Second coordinate
		vCoord[1] = coo;
		in >> coo;
		if (dim==3) vCoord[2] = coo;	// Third coordinate.
		
		// It is assumed that nodes do not exist previously.  Their names must be added at parentSystem
		parentSystem->addPointName(auxs,vCoord);
	}
	in >> auxs;		/** $EndNodes */
	in >> auxs;		/** $Elements */
	in >> nElements;
	connect.resize(nElements);

	for (i=0; i<nElements; ++i)
	{
		in >> auxi1;	// Element number (not used)
		in >> vertexTag;	// Element type (2=triangle, 4=tetrahedron)
		if (i==0)
		{
			if (vertexTag==2) nVertex = 3; else if(vertexTag==4) nVertex = 4; else cout << "**ERROR: Mesh non-consistent with dimension\n";
			vVertex.resize(nVertex);
			for (j=0; j<nElements; ++j) connect[j].resize(nVertex);
		}
		in >> auxi1;	// Integer not used
		in >> auxi1;	// Integer not used
		in >> auxi1;	// Integer not used
		in >> auxi1;	// Integer not used
		// Triangles
		in >> auxs;	// First vertex
		connect[i][0] = auxs;
		addSystemName(auxs,parentSystem->name);	/** System name is added to the point name */
		vVertex[0] = auxs;
		in >> auxs;	// Second vertex
		connect[i][1] = auxs;
		addSystemName(auxs,parentSystem->name);	/** System name is added to the point name */
		vVertex[1] = auxs;
		in >> auxs;	// Third vertex
		connect[i][2] = auxs;
		addSystemName(auxs,parentSystem->name);	/** System name is added to the point name */
		vVertex[2] = auxs;
		if (nVertex==4)
		{
			in >> auxs;	// Fourth vertex
			connect[i][3] = auxs;
			addSystemName(auxs,parentSystem->name);	/** System name is added to the point name */
			vVertex[3] = auxs;			
		}

		auxsv[0] = vVertex[0];
		auxsv[1] = vVertex[1];
		vComponents.push_back(new Truss(auxsv,area,matName,dim));
		
		auxsv[0] = vVertex[0];
		auxsv[1] = vVertex[2];
		vComponents.push_back(new Truss(auxsv,area,matName,dim));
		
		auxsv[0] = vVertex[1];
		auxsv[1] = vVertex[2];
		vComponents.push_back(new Truss(auxsv,area,matName,dim));
			
		if (nVertex==4)
		{
			auxsv[0] = vVertex[0];
			auxsv[1] = vVertex[3];
			vComponents.push_back(new Truss(auxsv,area,matName,dim));

			auxsv[0] = vVertex[1];
			auxsv[1] = vVertex[3];
			vComponents.push_back(new Truss(auxsv,area,matName,dim));
			
			auxsv[0] = vVertex[2];
			auxsv[1] = vVertex[3];
			vComponents.push_back(new Truss(auxsv,area,matName,dim));			
		}
	}
	in >> auxs;		/** $EndElements */
}

// ***
/*
void Lattice::writeFileGraph(ofstream &fGraph, string& gFormat)
{
		
}
*/

// ***

void Lattice::writeMSH(ofstream &fGraph)
{
	// edit to write only the surfaces	
	
	int i;
	int j;
	int nt;
	string auxs;
	int elementType = 0; 	// edit  number of the new element type for the gmsh file
	int nNewVertex = 0;  	// edit  number of vertexes of the new element 
	//char *outGraphFile;	// edit: get name of the graphic file  
	string numStep = "0"; 	// edit: current timestep 
	
	// Get the pointers to the points, to access the current coordinates
	vector<Point*>	vPoints(vPointNamesSurface.size()); // edit
	//vector<Point*>	vPoints(nNodes);	
	for(i=0;i<vPointNamesSurface.size();i++)
	{
		auxs = vPointNamesSurface[i]; // edit
		addSystemName(auxs, parentSystem->name);
		vPoints[i] = parentSystem->getPTPointLocal( auxs ); 
	}
	
	
	fGraph << "$MeshFormat\n";
	fGraph << "2 0 8\n";
	fGraph << "$EndMeshFormat\n";
	fGraph << "$Nodes\n";
	fGraph << vPoints.size() <<  "\n"; // edit
	for (i=0; i<vPoints.size(); ++i) //  edit
	{
		nt = vPoints[i]->globalTNum.size();	// Number of translational dofs
		fGraph << vPointNamesSurface[i]; // edit 
		for (j=0; j<nt; ++j)
		{
			fGraph << "\t" << parentSystem->getDof( vPoints[i]->globalTNum[j] ); 
		}
		if (nt<3) fGraph << "\t" << 0;	// Third coordinate, just in case is 2D
		fGraph << "\n";
	}
	fGraph << "$EndNodes\n";
	
	
	fGraph << "$Elements\n";
	fGraph << connectSurface.size() << "\n"; // edit
	
	if (vertexTag == 2) {elementType = 1; nNewVertex = 2;} // edit
	else if (vertexTag == 4) {elementType = 2; nNewVertex = 3;} // edit
	else {cout << vertexTag << " not (yet) available in lattice \n";} 
	
	for (i=0; i<connectSurface.size(); ++i)  // edit	
	{
		fGraph << i+1 << "\t" << elementType << "\t" << "3 " << "7 " << "6 "<< "0 ";  // edit	
		for (j=0; j<nNewVertex; ++j) // edit
		{
			fGraph << "\t" << connectSurface[i][j]; // edit
		}
		fGraph << "\n";
	}
	fGraph << "$EndElements\n";
	
	
	// get current time step from the outputfile name stored for that process
	// the graphic file name has to be: xxxxx-number.dat, number will be extracted from the file name, where "-" and "." are only used one time
	
	string outGraphFile(parentSystem->getOutGraphFile());			// get name of the graphic file
	outGraphFile = outGraphFile.substr(0,outGraphFile.find_last_of(".")); 	// get first part of the filename until ".(dat)"
	numStep = outGraphFile.substr(outGraphFile.find_last_of("-")+1);	// isolate current step
	
	
	fGraph << "$NodeData\n"; 
	
	fGraph << "1\n"; 			// one string follows
	fGraph << "\"Title or new name\"\n"; 	// name of the "view"
	
	fGraph << "1\n";			// one real tag follows
	fGraph << numStep << "\n";		// time value // can be different, should be changed to real time value !!!!
	
	fGraph << "3\n";			// three integer will follow
	fGraph << numStep << "\n";		// time step
	fGraph << "1\n";				// type: 1 = scalar field,...			 
	fGraph << vPointNamesSurface.size() <<"\n";	// number of associated nodal values
	
	// now start the nodal values
	for (int i = 0; i < vPointNamesSurface.size(); i++)
	{
		fGraph << vPointNamesSurface[i] << "\t" << "0.1" << "\n"; 	// "0.1" should be replaced by displacement or tension
	}
	
	fGraph << "$EndNodeData\n"; 
	
	/*
	// file with the gmsh command for the post processing animation
	// this is only a test a should be deleted later
	
	ofstream gmshcom;
	
	if (!(numStep.compare("0")))
	{
		gmshcom.open("gmsh_com.sh");
		gmshcom.close();
	}
	//double  test;
	//cin >> test;
	
	gmshcom.open("gmsh_com.sh",ios::app);
	
	if (!(numStep.compare("0")))
	{
		gmshcom << "gmsh -combine ";
	}
	
	string outGraphFilex(parentSystem->getOutGraphFile());	
	gmshcom <<  outGraphFilex << " -new ";
	
	
	gmshcom.close()	;
	
	// file with the gmsh command for the post processing animation
	// this is only a test a should be deleted later
	
	
	
	if (!(numStep.compare("0")))
	{
		gmshcom.open("gmsh_com_c.sh");
		gmshcom.close();
	}
	//double  test;
	//cin >> test;
	
	gmshcom.open("gmsh_com_c.sh",ios::app);
	
	if (!(numStep.compare("0")))
	{
		gmshcom << "sh post.sh" << "\n";
		gmshcom << "gmsh -combine ";
	}
	
	//string outGraphFilex(parentSystem->getOutGraphFile());	
	gmshcom << "c_" << outGraphFilex << " -new ";
	
	
	gmshcom.close()	;
	
	// file with the gmsh command for the post processing animation
	// this is only a test a should be deleted later
	
	ofstream post;
	
	if (!(numStep.compare("0")))
	{
		post.open("post.sh");
		post.close();
	}
	//double  test;
	//cin >> test;
	
	post.open("post.sh",ios::app);
	
	//if (!(numStep.compare("0")))
	 //{
	 //gmshcom << "gmsh -combine ";
	 //}
	
	//string outGraphFilex(parentSystem->getOutGraphFile());	
	post << "./post " << outGraphFilex << "\n";
	
	
	post.close()	;
	*/
}

// ***

void Lattice::writeVTU(ofstream &fGraph, bool& plotDispFlag, bool& plotProcFlag, bool& plotStressFlag)
{
	int i;
	int j;
	int nt;

	localFGraph << "<?xml version=\"1.0\"?>\n";
	localFGraph << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n";
	localFGraph << "	<UnstructuredGrid>\n";
	localFGraph << "		<Piece NumberOfPoints=\"" << vPPointSurface.size() << "\" NumberOfCells=\"" << connectSurface.size() << "\">\n";
	localFGraph << "			<Points>\n";
	localFGraph << "				<DataArray type=\"Float64\" NumberOfComponents=\"3\" format=\"ascii\">";
	for(i=0; i<vPPointSurface.size(); ++i)
	{
		localFGraph << "\n";
		localFGraph << "				";
		nt = vPPointSurface[i]->globalTNum.size();	// Number of translational dofs
		for (j=0; j<nt; ++j)
		{
			localFGraph << parentSystem->getDof( vPPointSurface[i]->globalTNum[j] ) << "\t"; 
		}
		if (nt<3) localFGraph << 0.0;	// Third coordinate, just in case is 2D
	}
	localFGraph << "\n";
	localFGraph << "				</DataArray>\n";
	localFGraph << "			</Points>\n";
	localFGraph<< "			<Cells>\n";
	localFGraph << "				<DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">";
	
	if(dim<3) nt=2; else nt=3;	// Hay que cambiarlo
	for (i=0; i<connectSurface.size(); ++i)
	{
		localFGraph << "\n";
		localFGraph << "					";
		for (j=0; j<nt; ++j) // edit
		{
			localFGraph << pNamePNumber[ connectSurface[i][j] ] << "\t";
		}
	}
	localFGraph << "\n";
	localFGraph << "				</DataArray>\n";
	localFGraph << "				<DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">";
	localFGraph << "\n";
	localFGraph << "					";
	for (i=0; i<connectSurface.size(); ++i)
	{
		localFGraph << nt*(1+i) << "\t";
	}
	localFGraph << "\n";
	localFGraph << "				</DataArray>\n";
	localFGraph << "				<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">";
	localFGraph << "\n";
	localFGraph << "					";
	if(dim<3) nt=3; else nt=5;	// Hay que cambiarlo
	for (i=0; i<connectSurface.size(); ++i)
	{
		localFGraph << nt << "\t";
	}
	localFGraph << "\n";
	localFGraph << "				</DataArray>\n";
	localFGraph<< "			</Cells>\n";
// Example of point Data
/*
	localFGraph << "			<PointData Scalars=\"scalars\">\n";
	localFGraph << "				<DataArray type=\"Float64\" Name=\"scalars\" format=\"ascii\">";
	localFGraph << "\n";
	for (i=0; i<vPPointSurface.size(); ++i)
	{
		localFGraph << 0.4 << "\n";
	}
	localFGraph << "				</DataArray>\n";
	localFGraph << "			</PointData>\n";
*/
	localFGraph << "		</Piece>\n";
	localFGraph << "	</UnstructuredGrid>\n";
	localFGraph << "</VTKFile>\n";		
}

// ***

void Lattice::writeOOGL(ofstream &fGraph)
{}

// ***

void Lattice::writeFileTemp(const double& time, ofstream & outF)
{}

// ***

void Lattice::setIniPosBody(vector<double>& iniValues)
{}